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ABSTRACT 

We study the dissipative effects of baryon physics on cosmic statistics at small 
scales using a cosmological simulation of a (50 Mpc/h) 3 volume of the universe. The 
MareNostrum simulation was performed using the AMR code RAMSES, and includes 
most of the physical ingredients which are part of the current theory of galaxy for- 
mation, such as metal-dependent cooling and UV heating, subgrid modelling of the 
ISM, star formation and supernova feedback. We re-ran the same initial conditions 
for a dark matter only universe, as a reference point for baryon-free cosmic statistics. 
In this paper, we present the measured small-scale amplification of a 1 and S3 due 
to baryonic physics and their interpretation in the framework of the halo model. As 
shown in recent studies, the effect of baryons on the matter power spectrum can be 
accounted for at scales k < 10 /i.Mpc - by modifying the halo concentration param- 
eter. We propose to extend this result by using a composite halo profile, which is a 
linear combination of a NFW profile for the dark matter component, and an expo- 
nential disk profile mimicking the baryonic component at the heart of the halo. This 
halo profile form is physically motivated, and depends on two parameters, the mass 
fraction fd of baryons in the disk, and the ratio Xd of the disk's characteristic scale 
to the halo's virial radius. We find this composite profile to reproduce both the small- 
scale variance and skewness boosts measured in the simulation up to k ~ 10 2 ft.Mpc -1 
for physically meaningful values of the parameters fd and A^. Although simulations 
like the one presented here usually suffer from various problems when compared to 
observations, our modified halo model could be used as a fitting model to improve the 
determination of cosmological parameters from weak lensing convergence spectra and 
skewness measurements. 

Key words: Gravitational lensing - cosmological parameters - Galaxy: disk - 
Hydrodynamics - large-scale structure of Universe 



1 INTRODUCTION 

One of the great challenges in modern cosmology is to under- 
stand the nature of dark energy, which is believed to domi- 
nate the energy budget in the universe (~ 70%) at low red- 
shift llRiess et al.ll 19981; iPerlmutter et al.lll999l ; ISpergel et all 



Astier et al.ll2006l ). Since the value of S1a directly im- 



120071 : 

pacts the rate of structure formation at recent epochs, the 
mass distribution and its time evolution bear the signature 
of the dark energy content of the universe. Cosmic shear 
measurements provide an independent method of probing 
the total mass distribution at large scales. Combined with 
photometric redshifts, it is even possible to extract the 3D 
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matter distribution and reconstruct the matter power spec- 
trum at different epochs. Comparing these measurements 
to theoretical predictions will set strong constraints on the 



cosmological parameters (e.g., 


Hu & Tegmarkll 19991: iHutererl 


2002: Amara & Refreeierll2006 


: Albrecht & Bernsteinll2007f). 



and among them both the equation of state w of dark energy 
and its possible evolution w' with redshift. 

The cosmic shear convergence power spectrum P K de- 
pends on the total matter power spectrum P, which con- 
tains the information about w and w' through the growth 
rate of structures. Extracting the dark energy equation of 
state from weak lensing signals therefore requires a good 
theoretical model for P, with a typical accuracy of a few 
per cent up to angular scales of about 10' (see f or exam- 
ple |Refregiei] [2003; [Bartelmamr^Schneidei] 200l|, for a re- 
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view). Substantial work has been done to measure and 
predict the dark matter power s pectrum from collision- 
less jV-body simulations (see, e.g., lEfstathiou fe Eastwood! 
Il98ll ; |jenkins et al.lll998h . Semi-analytic models for the dark 
matter power spectrum ha ve also been proposed, reaching 

et al.l 



the percent level accuracy llHamilton et al 



Il99ll :lja 

ll995l ; |Peacock fc Dodddll996l : ISmith et alj|2003t) . While the 
distribution of total matter is likely to closely follow dark 
matter at large scales, dissipative physics is expected to 
modify the power spectrum at small scales, and therefore 
possibly interfere with weak lensing measurements. 

The interest for the effects of baryons on the con- 
vergence power spectrum has led to the develop ment of 
semi- analytic halo models to estim ate effect of cold l|Whitd 
120041 ) and hot (jZhan fc KnoxH2004r i gas on the total matter 
power spectrum. More recently, numerical simulations have 
been carried out to complement th ose semi-analytical results 
l|jing et alJ200d : lRudd et al.ll200Sf ). While the exact effect of 
the baryons differ quantitatively between different models, 
the models and simulations agree qualitatively on a boost 
of the total matter power spectrum due to cold baryons at 
small scales. At k ~ 10 h.Mpc , this amplification has been 
found to be of around 10% at z — 0. Our theoretical under- 
standing of galaxy formation is however far from being com- 
plete. Current numerical simulations that include various 
complex baryons ph ysical processes suffer from the so- called 
overcooling problem l|Blanchard et al]|l992l ; [Cold [l991). with 
too many baryons condensing into gaseous and stellar disks 
when compared to observational constraints. Statistical ef- 
fects measured in Galaxy formation simulations, including 
the one used in the present work, are therefore likely over- 
estimated. If we can account for the effect of baryons at the 
required accuracy in this extreme case, real datasets will be 
probably even easier to deal with. 

We still need a flexible and accurate tool to account 
for the effect of baryons on the statistical properties of 
the matter density field in a parametrised model. The halo 
model has been developed in the last decade to meet these 
goals. The halo model is based on the idea that the mat- 
ter distribution in the universe can be described as a col- 
lection of individual halos, in which baryonic structures 
such a s galaxies form (|Nevman fc Scottlll952l ; IWhite fc Reesl 
1 1978ft . IScherrer fc Bertschingerl (|l99"lh ~ proposed a formalism 
to compute correlation functions of the continuous density 
field from a model of discrete virialized halos. Since then, 
there has been notable contribu tions and refinements to 
the ha lo model approach, such as iMa fc Frvl <l2000ri; ISeliakl 
(2000) and subsequent developments (see ICoorav fc Sheth 
|2002| for a review in the context of large scale structure) . 

As the halo model has proved to be a successful frame- 
work for describing statistical properties of the dark mat- 
ter density field in the non-linear regime, there has been 
also interest in extending it to baryons in the contex t of the 
Sunyaev-Zeldovich effect (Rcfrcaic r fc Tevssierll2002l ) and of 
the ga l axy d istrib ution (Scljak 2000). In previous studies, 
IWhitd l|2004 ) and IZhan fc Knoxl l|2004 ) have used the halo 
model with a baryonic component to describe the effect of 
cold and hot gas res pectively from a se mi-analytical stand- 
point. More recently. iRudd et al.l |2008r i have shown that the 
halo model can be used in a self-consistent way to describe 
the amplification of the power spectrum caused by baryons 
as measured in cosmological simulations. They proposed to 



modify the concentration parameter mass dependence of the 
dark matter halos to account for the collapse of baryons at 
small scale, leading to more concentrated halos. 

In this paper, we would like to extend the previous mod- 
els for cosmic statistics to smaller scales, where baryons are 
likely to dominate the total mass distribution. For that pur- 
pose, we use the results of a recent, high-resolution, cos- 
mological simulation, featuring state-of-the-art galaxy for- 
mation physics, thanks to the Horizon collaboration 1 . The 
simulation was performed on the MareNostrum supercom- 
puter at the Barcelona Su percomputer Centre using the 
RAMSES code l)Tevssierll2002l l. including a detailed treatment 
of metal-dependent gas cooling, UV heating, star formation, 
supernovae feedback and metal enrichment. 

The simulation data are compared to the analytical pre- 
diction of a modified halo model, taking into account small 
scale baryons physics in an ad-hoc way by adding to the 
halo mass profile a small baryonic component, modelled as 
an exponential disk with mass fraction and scale length as 
the only 2 additional free parameters. This approach, which 
modifies the shape of the halo profile, is in essence similar 
to the one of lWhitel (|2004l 'l. which we use as a starting point 
for our theoretical model to compare against our numerical 
simulation. In contrast to previous studies, we also compute 
the effect of baryons on the skewness of the mass distribu- 
tion. It has been shown that measuring the third moment 
of the cosmic shear is of paramount importance, since it can 
break the degeneracy in the cosmological parameters esti- 
mation based on the power spectrum alo ne, and reduce the 
corresponding error bars by a factor of 2 dBernardeau et al.l 
1 19971 ; I Jain fc Van Waerbekd I2OO0I ; iTakada et al.ll2000l l. Us- 
ing only the two additional parameters of the model, we 
can fit the simulation data with great accuracy, for both 
the power spectrum and the skewness. This has important 
consequences for future weak lensing surveys, since the disk 
parameters of the model could be fitted together with the 
cosmological parameters, promoting baryons physics from 
a mere systematic effect to an additional probe of the un- 
derlying cosmolo gical model. With in the modified concen- 
tration model of lRudd et al.l (120081) , statis tical bias effects 
have been studied by Zentner et alj (|2008T l , and further by 
iHearin fc ZentaeH (2009) in the context of the test of general 
relativity by weak lensing surveys. 



2 STATISTICAL ANALYSIS OF THE 
MARENOSTRUM SIMULATION 

2.1 The MareNostrum simulation 

We have performed a cosmological simulation of unprece- 
dented scale, using 2048 processors of the MareNostrum 
computer installed at the Barcelona Supercomputing Centre 
in Spain. We h ave used intensively the AMR code RAMSES 
l|Tevssierl [2002) for 4 weeks dispatched over one full year. 
This effort is part of a consortium between the Horizon 
project 1 in France and the MareNostrum galaxy formation 
project in Spain 2 . The originality of this project relies on 
using a lot of (if not all) physical ingredients that are part 
of the current theory of galaxy formation, and at the same 
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time cover a large enough volume to provide a fair sample 
of the universe, especially at redshifts above 1. 

The physical processes we have included in our simula- 
tion are described now in more detail. We have considered 
metal-dependent cooling and UV hea ting using the Hard t 
and Madau background model (see lOcvirk et al.l [2008). 
We have incorporated a simple model of supernova feed- 
back and metal enrichment using t he implementation de- 
scribed in iDubois &: Tevssierl l|2008T ). For high-density re- 
gions, we have considered a polytropic equation of state 
to model the co mplex, multi-phase and turbulent struc- 
ture of the ISM |Yepes et al.l Il997l; ISpringel fc Hernqui; 



20031') in a simplified form (see lSchave fc Dalla Vecchiall200S 
Dubois fc Tevss icr 2008): the ISM is defined as a gas with 
a density above no ~ 0.1 H/cc. Star formation has also 
been included, for ISM gas only (uh > no), by spawn- 
ing star particles at a rate consistent with the Kennicutt 
law derived from local observations of star-forming galaxies. 
In more mathematical terms, we have p\ = p gas /f* where 
i* = (tih /no) -1 to and to = 8 Gyr. Recast in units of the 
local free-fall time, this corresponds to a star formation ef- 
ficiency of 5%. The simulation was started with a base grid 
of 1024 3 cells and the same number of dark matter parti- 
cles, and the grid was progressively refined, on a cell-by-cell 
basis, when the local number of particles exceeded 10. A 
similar criterion was used for the gas, implementing what 
is called a quasi-Lagrangian refinement strategy. Five ad- 
ditional levels of refinement were considered, providing a 
resolution between 1 and 2 kpc physical at all times. In this 
way, our spatial resolution is consistent with the angular res- 
olution used to derive the Kennicutt law from observations. 
On the other hand, we are not in a position to resolve the 
scale height of thin cold disks so that the detailed galactic 
dynamics might be affected by resolution effects. 

The simulation was run for a ACDM universe with 
fi m = 0.3, n A = 0.7, Q b = 0.045, H = 70 km/s/Mpc, 
as = 0.9 in a periodic box of 50 Mpc/ft. Our dark matter 
particle mass (m p ~ 8 x 10 6 Mq /h) , our spatial resolution (1 
kpc physical) and our box size make this simulation ideally 
suited to study the formation of galaxies within dark mat- 
ter halos, from dwarf- to Milky- Way-sized objects at high 
redshift. For large galaxies, we can nicely resolve the radial 
extent of the disk (not its vertical extent), while for small 
galaxies, we can resolve the gravitational contraction of the 
cooling gas, but barely the final disk. The simulation was 
stopped at redshift z ~ 1.5 because we ran out of allocated 
time. The total number of star particles at the end of the 
simulation was above 2 x 10 s , and the total number of AMR 
cells was above 5 x 10 9 . 

To quantify the effects of baryons on statistical proper- 
ties of the mass distribution, the MareNostrum run, which 
includes dissipative physics, was re-run from the same ini- 
tial conditions with baryon mass converted to dark matter. 
This latter dark matter only (henceforth DMO) simulation 
serves as a reference for statistical quantities without the 
presence of dissipative physics. Both runs were carried out 
up to redshift 2 = 2, which we will consider in the rest of 
this paper. It is already late enough to witness interesting 
structures such as well-formed galaxy disks. 
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2.2 One-point statistics 

Meaningful statistics of the density field can be extracted 
from different statistical quantities, such as the n-point cor- 
relation functions, the density PDF, or the one-point cumu- 
lants. By far, the easiest quantities to measure are the one- 
point moments S P (R), i.e. the p-ih order cumulant of the 
smoothed density field as a function of the smoothing scale 
R. The S p parameters have also been studied extensively, 
whether from a theoretical st andpoint |Balian fc Schaefferl 
1989; Szapudi & Szalay 1993), in the perturbation theory 
framework ( Bernard eau 1994), or in simulations and obser - 
vations (see, e.g., Colombi et al. . 2000; Marinoni et al . 2008). 
For the 50 Mpc/h box of MareNostrum, we have computed 
the statistics for scales ranging from 15 kpc/h to 2 Mpc/ft. 

With weak lensing applications in mind, we are primar- 
ily interested in the total mass statistics. In the case of the 
dissipative simulation, this requires a consistent treatment 
of both dark matter particles and gas cells. 

The total local density in the dissipative simulation can 
be written 



Ptot = pg + Pdm + p s , 



(1) 



where p g , Pdm and p s are the local gas, dark matter and star 
densities respectively. However, because of the different na- 
ture of the gas (which is a continuous cell-based quantity), 
and stars and CDM (which are modelled as collisionless par- 
ticles), the three matter components should be treated sepa- 
rately. One could simply evaluate pnM and p a by binning the 
particles into cells to obtain a local estimate of the densities, 
and then simply calculate p t ot as in Eq. [T]and computing its 
moments. However, as we discuss below, the discrete nature 
of particles mandates a special treatment, and we chose in- 
stead to compute the moments and cross-correlations of the 
different species separately, and then reconstruct the mo- 
ments of the total density field as we now describe. 

Obtaining the moments of the gas distribution involves 
no theoretical difficulty. The gas density of the whole sim- 
ulation box is mapped onto a 2048 3 grid from the AMR 
cells using a donnercell scheme, where the resulting value in 
each cartesian cell is directly copied from the finest AMR 
cell covering it. To determine the moments of the smoothed 
gas density field for a given comoving smoothing radius R, 
we compute the average of the density over cubic regions 
of volume |7rJ? 3 . We restrict ourselves to values of R corre- 
sponding to smoothing boxes which span an integer number 
of base grid cells. The average densities in such cubic re- 
gio ns are computed usi ng a fast convolution algorithm (see 
e.g. iBlaizot et al.ll2"006h . and the moments over all such re- 
gions are then evaluated. Since the simulation directly pro- 
vides the continuous gas density p g , this prescription yields 
unbiased estimates of the moments of the gas distribution. 

Particles require a somewhat more careful treat- 
ment. The statistics of particle distributions are readily 
studied using a counts-in-cells analysis (see for example 
Balian fc Schaefferlll989l: iBouchet fc Hernquist|[l992l ; IShethl 



19961) ■ The idea is to count particles within the same cu- 
bic cells of scale R used for the smoothing of the gas den- 
sity. The particle counting is implemented by first bin- 
ning the particles in to the base grid using a ne arest grid 
point (NGP) scheme (|Hocknev fc Eastwood|[T98lh . and then 
counting particles in the cubic regions, again by using fast 
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convolution. This is indeed equivalent to computing a local 
particle density by NGP, and then performing the _R-scale 
smoothing. In this case, however, simply computing the mo- 
ments of the resulting data will introduce shot n oise effects 
|Szapudi fc Szalay|[l993l ; iBernardeau et alj|2002h . 

It is possible to correct for these effects using factorial 
moments. Let us consider a continuous field p sampled by a 
finite collection of particles. Given a cell of size R and volume 
v — R s , we call p — i J p(x)d 3 x the average value of p 
over the cell. Equivalently, p is the value of p smoothed at 
scale R at the centr e of the cell. The local P oisson sampling 
hypothesis (see e.g. IBernardeau et aTl l2002) states that the 
distribution of the number N of particles in the cell follows 
a Poisson probability mass function of mean pv / m p , where 
m p is the mass of a single particle. Letting 

(N) h =N(N-l).--(N-k+l), (2) 

the factorial moments are defined as 

Fk = ((N) k ) p = (N(N - 1) • • ■ (N - k + l)) p , (3) 

where N is the cell particle count, and (. . .) p denotes the 
average over the Poisson distribution of the part icle sam- 
pling. It has been shown (|Szapudi fc Szalavlll993t ) that the 
Ffc yield unbiased estimators of the moments of the under- 
lying smoothed density field p at the scale of the cell size, 
in the sense that 

Let us now consider the density fields smoothed at scale 
R for the gas, dark matter and stars, p g , pdm and p s re- 
spectively. For the sake of readability, we shall drop the 
tilde notation, and the density fields are to be understood as 
smoothed at the scale R in the rest of this section. Using Eq. 
[T] we can express the moments of the smoothed total den- 
sity field ptot as a function of the moments and correlations 
of the individual species using the multinomial theorem: 



We can finally compute the Sk parameters, which are 
defined as 



-fc 

P = 



k 

ki, fe, 



fc 



(p^PdWv 3 )- (5) 



In this equation, {• • •) denotes the ensemble average over 
all realisations of the underlying density fields, not to be 
confused with the average (•■•)□ over particle samplings for 
a fixed realisation of p introduced in equation O 

Provided we can compute all correlations of the form 
{pg 1 /°dm/°s 3 ) ' we are now m position to reconstruct the mo- 
ments of the smoothed total density field. Under the local 
Poisson sampling hypothesis, one can deduce from Eq.[3]the 
identity: 

/ fci k 2 k 3 \ /m DM \ fc2 / m s \ fc 3 , k . x 
<P/Pd>, 3 ) = [— ) { — ) W (Nnu) k2 (N s ) k3 ) , 

(6) 

which involves the definition ([2]), and where 771dm and m s 
are the dark matter and star particle masses. Since the Pois- 
son processes of the counts-in-cells for the different particle 
species are independent of each other, Eq. [5] involves no ap- 
proximation, even though the underlying fields pdm and p 3 
are correlated. 

From the moments ([5}, we can straightforwardly com- 
pute the moments of the total matter overdensity (S^^ = 

((ptot/ptot — from the binomial theorem. 



S k (R) 



(7) 



where the c subscripts denote the connected moments of 
the smoothed density field, whose generating function is the 
logarithm of the generating function of the (<5 fe ). 

2.3 Power spectrum 

Because of the particular significance of the 3D total matter 
power spectrum P(k) in the convergence power spectrum, 
we have also measured P(k) in the dissipative and DMO 
simulations, in addition to the one-point statistics. The vari- 
ance of the total matter density field smoothed at scale R 
can be expressed as: 

a2{R) = 2^J X fc3p(fe) l^ fei? )| 2 > ( § ) 

where W is the Fourier transform of a spherical top-hat 
window function with volume unity: 

3 , . 



W(x) 



(9) 



Various sophisticated techniques for estimating the power 
spectrum have been proposed, es pecially for correcting mass 
assignment and sam pling effects l|jindl200^ : ICui et all 2008; 
IColombi et al. 2009). Since the 2-point correlation function 
(or, equivalently, the power spectrum) is not our primary 
interest in this paper, we have settled for a simple method 
which we expect to yield reasonable results, even if not as 
accurate as our one-point moments measurements. 

The gas and particles densities were mapped onto a 
2048 3 base grid and added up, using donnercell for the 
gas and NGP binning for the dark matter particles. The 
spectrum is computed usin g FFT folding (see I Jenkins et ail 
ll998l : IColombi et al-lfiopgh and corrected for the NGP con- 
volution and shot noise bias effects (|Hocknev fc Eastwood! 



2.4 Results 

Because of cooling, the baryons will condense to form dense 
structures such as disks at the centre of dark matter ha- 
los. Figure [T] shows a density map of one of the biggest 
MareNostrum halos, together with contours of the density 
ratio Pbar/pcDM- The effect of cooling can be seen as dense 
baryon-dominated regions at the cores of the halos and halo 
substructures. 

The small-scale baryonic features directly impact the 
density statistics at small scales: as the smoothing scale de- 
creases, the disks become more and more apparent in the 
density PDF as peaks in the high-density regions. We can 
expect this process to broaden the distribution, thereby in- 
creasing the variance, and as only the higher-density regions 
are affected, the skewness should also increase. 

The computed variance a 2 and skewness S3 from the 
MareNostrum dissipative and DMO simulations is presented 
on Fig. [2] Comparing the DMO simulation (solid black) with 
the total matter in the dissipative run (blue dashes), we in- 
deed note that the presence of baryonic physics dramatically 




Figure 2. Variance (left) and skewness (right) of the smoothed density field of different species at z = 2, as a function of the smoothing 
scale in the MareNostrum dissipative and DMO simulations. The solid line shows tr 2 and S3 for dark matter in the DMO simulation, 
while the dashed and dotted lines correspond to the dissipative simulation: short dashes for dark matter, long dashes for total matter, 
and dots for baryons. The error bars on the DMO data are one-sigma wide and determined by the subvolumes method as described in 
the text. 



amplifies both a 2 and S3 at small R. At k ~ 10 Ti.Mpc -1 , 
the power spectrum boost reaches about 35% (see Fig. Q, 
most of which is caused by cold baryons (stars). Because 
our study is carrie d out at z = 2, p r ecise comparisons with 
previous results of Ijing et ail (|200ot ); iRudd et al.l (120081 ) are 
difficult. Note however that we observe the same qualita- 
tive effects. The variance plot on Fig. [5] also demonstrates 
the p resence of CDM adiabatic contraction (|Gnedin et al.l 
|2004| ). As the gas cools down, its contraction drags the dark 
matter into local potential wells created by dense baryon 
clumps. This effect results in a net condensation of the dark 
matter, whose effect on variance can be seen by comparing 
the DMO run (solid black curve) with CDM of the dissi- 
pative run (short-dashed curve). Both observed boosts and 
dark matter contraction effec ts are well in acc ordance with 
the results presented in Weinberg et all |2008l ). 

Because of the relatively small size of the MareNostrum 
simulation box, the results presented on Fig. [2] are contam- 
inated to some degree by cosmic variance and finite volume 
effects. We have estimated those effects in the MareNostrum 
DMO simulation. Note that the rigorous determination of 
error bars is beyond the scope of this article, and we do not 
expect baryons to modify those uncertainties significantly. 

The cosmic variance and finite volume effects on the 
statistical quantities were sampled by three different inde- 
pendent methods. We have run a set of 10 smaller 256 3 
cosmological simulations up to z = 2 with the same box size 
and power spectrum as the MareNostrum box, only with dif- 
fering random phases. The statistical quantities were then 
computed on each box, and the variance of those quanti- 
ties over the 10 boxes were used as a first estimate of the 
MareNostrum cosmic variance effects. While such ensemble 



simulations are easy to carry out, this method is known to 
underestimate the actual cosmic variance, as all the reali- 
sations of the initial density field are constrained: the total 
box matter density is fixed to the background matter den- 
sity of the universe. In addition, this method cannot be used 
to determine the variance at small scales because of the low 
resolution of the ensemble simulations. Relative cosmic er- 
ror derived from this set of simulations is presented on Fig. 
[3] (das hed curve). The FORCE code (FORtran for Cosmic 
Errors Colombi fc _Szapudi|[200ll ), implements the results of 
ISzapudi et all (| 19991 ) and provides cosmic variance estima- 
tions given the values of the density cumulants. The corre- 
sponding cosmic error, based on the MareNostrum DMO cu- 
mulants, is shown as the solid curves on Fig. [3] This estima- 
tion relies on a perturbative expansion which breaks down 
when relative errors approach unity. As the MareNostrum 
errors range from about 5% to 30%, the FORCE computa- 
tion still holds, but the quality of the estimation is impacted, 
especially at small scales where the errors on high-order cu- 
mulants increase. To confirm the FORCE results at small 
scales, we have studied the variance of the statistical quan- 
tities over a r andom samp le of cubic subvolumes of size i. Let 
€x{£,R) = \J var(X(R))/A"(_R) be the relative cosmic error 
of a statistical quantity X at scale R defined on subvolumes 
of size I. To obtain the cosmic variance of the whole simula- 
tion box (i.e., ex(L,R) for all R), we computed ex(£,R) for 
i — L/8, L/16, L/32 and extrapolated in £ to £ = L assum- 
ing the power-law form txil, R) = £x(L, R)(£/L) v . This last 
estimation of the error is represented on Figure [3] in dotted 
lines. None of these methods ensures accurate determination 
of the errors over the whole range of scales, however, they 
paint a clear picture of cosmic variance in the MareNostrum 
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500 kpc//i 



Figure 1. Map of the projected dark matter density of one 
of the largest halos in the MareNostrum simulation (iJ v ir = 
0.59 Mpc//i comoving shown as the dashed circle) at z = 2. 
The contours represent isovalues of the baryon to dark matter 
density ratio Pbar/PCDM- The outer black contours correspond to 
Pbar = O.lpcDMi while the inner red contours delimit equal den- 
sities regions. The total matter density is baryon-dominated at 
small scales well within the halo core. The bright central galaxy 
clearly stands out of the halo substructures, whose distribution 
within the halo remains mainly unaffected by the presence of 
baryons ( see IWeinberg et alj|2008h . 



simulations. As can be seen on Fig.0 the observed boosts in 
a 2 and S3 are well above cosmic variance effects. Note that 
scales comparable to the MareNostrum box size correspond 
to a patch of z = 0.5 sky extending over about 4 squared 
degrees. 

For our present study, it is important to notice, how- 
ever, that since both the DMO and dissipative runs have 
been performed using the same set of random phases for the 
initial conditions, they suffer from the same such effects. As 
a consequence, the corresponding errors in the two runs are 
strongly correlated, and ratios of statistical quantities such 
as ""tot/o'DMO are mostly devoid of finite volume contam- 
ination. For the rest of this paper, we will therefore only 
consider such amplification ratios (or "boosts") for the vari- 
ance and skewness of the total matter in the dissipative run 
with respect to the dark matter of the DMO run. 



3 A MODIFIED HALO MODEL FOR 
BARYONS 

3.1 The halo model 

The halo model provides a well-tested and flexible frame- 
work for the study of the properties of matter distribu- 
tion in non-linear stages of gravitational collapse. While 
first studied in the c ontext of the galaxy distribution 
|Nevman fc Scottl Il952l ) , it has become a full-fledged and 



now mature tool for cosmological statistics through signifi- 
cant contributions and improvements to its various ingredi- 
ents. 

Attempting to reproduce the effect of baryons on 
the boost factors of variance and skewness requires us to 
model both the D MO and dissipative matter distributions. 
iRudd et al.l (|2008l ) have shown that modifying the halo con- 
centration relation can account for the variance amplifica- 
tion at scales k < 20 h.Mpc~ . In this paper, we will use a 
standard halo model to describe the dark matter of the DMO 
run. We base our halo profile for the total mass on the DMO 
halo model, but instead of modifying only c(M), we propose 
to also modify the halo profile itself. As discussed previously, 
the quantity of interest is the boost of the statistics (i.e. the 
amplification of the variance and skewness witnessed on the 
total matter halo model with respect to the reference halo 
model). We now briefly describe the different key ingredients 
which take part in the computation of cr 2 {R) and Ss(R) in 
the halo model. 

Statistical description of the density field as a set of 
virialized halos requires the specification of the mass distri- 
bution of the halos (the mass function), their density profiles 
and associated mass parametrisation, and a model for halo- 
halo correlations. 

A simple model f o r the halo mass function was given 
by IPress fc Schechterl (|l974l ) based on the spherical col- 
lapse model. Since then, there has been more convincing 
derivations of the Press-Schechter result, as well as attempts 
to ta ke into account non-symmetric collapses and tidal ef- 



fects |Bond et all 



ISheth fc Tormen 



1991 



2002) 



lAudit et al.lll997l ; ISheth et alJboOll ; 

resulted in other 



These studies 

>arametrizations, such as the Sheth- Tormen mass function 
Sheth fc Tormenlll999h . 

In this study, we use the Sheth- Tormen prescription for 
the halo mass function, as it turns out to be a better fit to 
our simulations than the Press-Schechter form. In normal- 
ized units, the Sheth- Tormen mass function reads: 



n{m)dm = A(p)J^ (l + (qv 2 ) P ) 



(10) 



x i^exp 



where v = 5 c /a. S c « 1.68 is the collapse density threshold in 
the spherical collapse model, and p « 0.3, A(p) « 0.322, q ~ 
0.75. 

We have also introduced a mass cutoff in the halo mass 
function to account for the small box size of the MareNos- 
trum simulation. The value of the cutoff is chosen slightly 
above the mass of the biggest halo found in the simulation, 
which is around 5.10 13 Mq. While the cutoff has little effect 
on the variance as computed by the halo model, the skew- 
ness drops significantly at large scales, resulting in a bet- 
ter fit against the measured 53. This is not surprising since 
high-order moments at larg e scale are sensitive to rare events 
(such as massive halos, e.g. IColombi et al] [T994). which are 
not well sampled by the MareNostrum box. 

For the DMO model halo p rofile, we use the standard 
NFW form (jNavarro et al.lll997r ): 

_ c(M)r 



unfw (r|Af) ccx 1 (1 + a;) 



Ry 



(11) 



and u(r\M) is normalized such that J u(r\M) d 3 r = 1. Our 




Figure 3. Estimates for the relative cosmic errors Aa 2 /a 2 and AS3/S3 for each method described in the text. The dashed curves 
correspond to the 10 ensemble simulations, the solid curves to results of the FORCE code, and the dotted curves to the subvolumes 
estimation. 



halo virial radius i? V i r is defined such that, for a halo of mass 



M, 



we have M = ^npR* 



.A, with A = 200. Note that the 



mass M of a halo is related to its comoving Lagrangian size 
R by M = 1-irpR 3 . The NFW model has proved to fit nu- 
merical dark matter profiles over a large range of masses and 
radii with some dispersion i n the conc entration parameter 
c(M) l|Kravtsov et al.lll99g| ; Ijind l200Ch . The central loga- 
rithmic slope of dark matter profiles , which is —1 in the 
case o f NFW, is currently debated fsee lFukushige fc Makinc 
1997; Moore ct al. 1998, and more recently ISpringel et al. 
2008; Stadcl ct al. 2009). Note however that in the presence 



of dissipative physics, baryons are likely to affect the inner 
slope. 

The concentration parameter c (M) is parame t rized in 
our model according to the result of iBullock et a l. (1999): 



c(M, 



C() 



CO 



9, b', 



-0.13. 



(12) 



1 + Z 

This power-law model has been found to be a good fit 
to numerical simu lations also for dark energy cosmologies 

i|Dolag et al.ll2004Tl. 

Following Scherrer fc Bertschingerl (|l99ll ). we can ex- 
press the density 2-point correlation function £(r) as: 



£(r)=£ifc(r)+&*(r), 



(13) 



where £1^ represents the contribution to the correlation 
function from mass within the same halo, and f^fa contains 
the contribution from mass located in different halos. 

£ih is essentially the autocorrelation of the halo profile, 
and its contribution to the power spectrum is: 

Pih{k) = J n{m) hj\ \u(k\m)\ 2 dm, (14) 

where u(k\m) is the Fourier transform of the halo profile for 



a halo of mass m. We comp ute the halo-halo contribution 
followi ng Mo & White ( 1996) and its sub s equent extensions 
ilMo et al.1 Il997l ; ISheth fc Lemsonl 1 19991 : ISheth fc Tormenl 
Il999h . Assuming deterministic biasing on large scales, we 
can write the £hh contribution from two halos of masses Mi 
and M2 as: 



i hh {r\Mi,M2) 



b(Mi)b{M 2 )£(r) 
6(Mi)6(M 2 )£ lm (r) 



(15) 



where £i in is the matter correlation function from linear the- 
ory. This prescription is accurate at large scales, and consis- 
tent with the choice of mass function pr ovided the bias b (M) 
is computed from f(u) as prescribed in I Mo et al.l (| 19971 ). 

Now in possession of a halo model for £(r) (and there- 
fore its Fourier transform P(k)), we can evaluate cr 2 (R) us- 
ing Eq. U 



3.2 Skewness in the halo model 

While in principle the halo model ingredients presented so 
far fully determine the statistics of the density field, addi- 
tional work is needed to extract Ss(R). 

At large enough scales, the one-point st atistics Sk 
may be computed using perturbation t heory (fFrvl Il984l: 
Juszk iewicz et a l. 1993; Bernardeaull 19941 : iBernardeau et al.l 
2002TI . which yields 



qPT 34 

S3 = — +-y, 



(16) 



where 7 = din a 2 (R)/d In R. However, in the MareNostrum 
simulation at z = 2, PT is only expected to be valid at scales 
greater than ~ 1 Mpc//i. A first interesting refinement tak- 
ing discrete halos into account is the Poisson cluster model, 
where halo-halo correlations are neglected and profiles are 
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assumed to be point-like (|Shetblll996h . Halo profiles, how- 
ever, are responsible for most of the behaviour of small- 
scales statistics, and thus neither perturbation theory and 
the point-cluster model are appropriate for our study. 

Fortunately, the full computation of the higher- 
order cumulants Sk in the halo model was developed in 
IScoccimarro et all ([2001). Following the authors, we define: 



u m (R,v) 



k 2 Ak 
2tt 2 



u(k\v)) m \W(kR)f 



(17) 



3+1 



(18) 



where R is such that 5 c /<j(R) = v. In these notations, the 
third cumulant of the density field in the halo model writes 



(S 3 ) c = Sf T a? in + 3af in Aifi + A 0> 



(19) 



3.3 Halo model results 



We have tested some families of halo profiles to attempt to 
reproduce the observed effect of baryons on the statistics 
of the density field. The reference halo model for the DMO 
simulation is based on a NFW profile with t he commonly 
used c(M) relationship of lBullock et all l|l999f ) as written in 

Eq.na 

As suggested by previous numerical studies |Rudd et all 

2008), an increase in Co and a steeper concentration slope 
6 are expected to reproduce - at least partially - the in- 
crease in power at small scales due to baryonic physics and 
radiative processes in particular. We have accordingly tried 
to adjust the concentration parameters with a NFW profile 
to obtain a good match for both the variance and skewness 
at small scales. The power spectrum, variance and skewness 
bo osts for a NFW-ba sed model with parameters comparable 
to lRudd et ail (|2008h (c = 20, & = -0.15) are presented as 
the dotted curves on Fig. [4] and [5] This model reproduces 
the MareNostrum variance and power spectrum boosts down 
to a scale of about 0.5 Mpc/h. At smaller scales however, 
the halo model underestimates the variance amplification. A 
large part of this discrepancy is likely due to the difference 
in the simulation codes and physical modelling between the 
two studies. Note however that the skewness S3 of this halo 
model lacks much of the measured small-scale amplification, 
as can be seen on Fig. [S] The distinctive bend is also not re- 
produced at all, which suggests the profile form distributes 
matter too evenly across scales. 

With the partial success of this profile, one might ex- 
pect NFW profiles with higher concentrations to yield better 
fits. It turns out however that reasonably fitting the vari- 
ance boost at small scale requires very high values of Co, 
exceeding 30. Such high values of the concentration param- 
eter are too high to be accepted as physically meaningful. 
Yet more importantly, while increasing Co will indeed boost 
the variance, it fails to reproduce at all the corresponding 
small-scale skewness amplification. This can be seen on Fig. 
[S] and the S3 boost of a pure NFW halo model remains 
essentially flat for varying values of cq, with a very weak 
dependence. 

This leads us to believe that, while the NFW pro- 
file with adjusted concentration parameters has merits in 




10- 1 



10° 10 1 
k [/iMpc -1 ] 



10 2 



Figure 4. Relative power spectrum amplification due to baryons 
at z = 2. The solid curve is the measured power spectrum, the 
dotted curve is a NFW profile with cq = 20, b = —0.15, and the 
dashed curve is the halo model with our composite halo profile. 




Figure 6. Error on the amplification of the 3D total mass power 
spectrum at 2 = 2 for the halo models represented as a function 
of the angular mode i in the flat sky approximation. The dashed 
curve is the reference DMO model (i.e. without any boost), the 
dotted curve is the pure NFW model with modified c(M), and the 
solid curve is the composite halo model amplification. The light 
and dark shaded areas are estimates of the expected experimental 
errors on C'i for the CFHT Wide-field and LSST experiments, 
respectively. 
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R [Mpc//i] R [Mpc/h] 

Figure 5. Effect of baryons on the variance and the skewness 53 boost factors, as measured on the MareNostrum simulation (solid 
curve) and modelled by a NFW profile with cq = 20, b = —0.15 (dotted curve), and the composite profile (dashed curve). 



modelling the variance amplification caused by dissipative 
physics, it can only paint a limited picture of the statistical 
properties of the density field in the presence of baryons. 
As increasing Co essentially amounts to concentrating more 
matter within the central region of the halos, we naturally 
turn to other centrally-concentrated halo profiles. 

One way to concentrate more matter within the centre 
region is by using families of profiles with steeper central 
cusps than NFW of the form: 

u a (x) =x' a (l+x) a - 3 , (20) 

where a = 1 yields an NFW profile. We have tested this 
family of profiles on a wide range of values 1 < a < 2.5. For 
each value of a, we attempted to find an best-fitting value 
of (co, b), again by exploring the parameter space. It is inter- 
esting to note that high values of a, in the range [2.0, 2.15], 
produce to some extent both the a 2 small-scale steepen- 
ing and a strong 5*3 amplification. Isothermal (a = 2) pro- 
files are known to be a good description of th e total density 
in haloes hosting elliptical g alaxies (see e.g. iGavazzi et al.l 
120071 : Ivan de Ven et all 12009 ). In the case of the MareNos- 
trum simulation however, this property seems coincidental, 
as the simulated physics form no truly elliptical galaxy com- 
parable to observations. Moreover, the residuals of the best 
a 2 and S3 fits for such profiles cast doubt on the legitimacy 
of the analytical form u a for the statistical analysis of the 
simulation. 



3.4 A modified halo profile 

A good candidate profile which is both centrally- 
concentrat ed and physically-motivated is a composite halo 
profile (see lWhitell2004l ; IZhan fc Knoxll2004) . parametrized 



by the dimensionless parameters fd and Ad: 

«/d.Ad(r|A0 = (1 - fd) UNFw(r|M) 
+ fdU C x P ,\ d (r\M) , 

where u C x P ,A d is a spherically averaged exponential disk pro- 
file with length scale r<j proportional to the halo's virial ra- 
dius: 

« cxp ,a>|M) cx eXp( ~ r/r<i) , r d = X d R vlT . (22) 

The dimensionless parameter Ad is essentially the spin pa- 
rameter of the halo, and defines the disk scale r<j. The profile 
Uf dt \ d features a central r _1 cusp and behaves like the NFW 
profile for radii bigger than the disk length scale rd- How- 
ever, because of the profile normalization, it concentrates 
more mass within the central exponential than a pure NFW. 
Uf dl \ d can be seen as a halo profile concentrating a fraction 
fd of the mass within a central exponential disk profile, and 
the remaining 1 — fd in a standard NFW component. 

This form of composite profile is physically motivated. 
The total mass distribution in group-sized halos is known to 
be well described by a halo component and a concentrated 
comp onent correspo nding to the bright central galaxy (see, 
e.g.. lDubin ski 1998). The presence of baryons does not fun- 
damentally change the diffuse halo component: the distri- 
bution of satellite galaxies within halos is very similar to 
the halo occupation distribution of dark matter substruc - 
tures in pure TV-body simulations (see Wei nberg et al . 2008). 
This suggests keeping a NFW profile to account for the dark 
matter, diffuse gas and halo substructures, while introduc- 
ing a spiked central component mimicking the bright central 
galaxy's disk. We may expect this NFW profile to be more 
concentrated than in the dark matter only case, because 
of the adiabatic contraction of the CDM due to the pres- 
ence of baryons. This will therefore lead to an increase of 
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Co in the c(M) relationship of equation (|12p . For the com- 
posite profile, fd is to be understood as the fraction of the 
total halo mass which resides in the galactic disk in the 
form of baryons. As most formed galaxies found at z — 2 in 
MareNostrum simulation are spirals, we restrict ourselves in 
this paper to an exponential disk profile for the central com- 
ponent. We believe this form captures the essential features 
of the dense central baryonic regions which are important 
for the halo model. It also places interesting constraints on 
the profile parameters fd and Xd, as mass fractions and an- 
gular momenta of disks are well-studied, both theoretically 
and observationally. We further assume both Ad and fd to be 
independent of halo mass. The assumption that the disk size 
is a fixed fraction of 7? v i r corresponds to the singular iso ther- 
mal sphere model of disk formation fsee lMo et al.lll998h . We 
postpone refinements of this model to future work. 

Here again, we explored the (fd, Ad) parameter space to 
find a good fit to the MareNostrum data. Our best model 
has parameters: 

{c = 13.5 
b = -0.15 
fd = 0.09 (26) 
A d = 0.025 

The corresponding power spectrum, variance and skewness 
boosts are represented on Fig. U and [5] as dashed curves. 
This halo profile reproduces accurately both the measured 
a 2 and S3 amplifications down to the smallest scales. 



4 DISCUSSION AND CONCLUSION 

With a base grid resolution and particle count of 1024 3 and a 
box size of 50 Mpc/ft, the MareNostrum simulation resolves 
the length and mass scales of galactic disks while also pro- 
viding a volume large enough for cosmological studies. This 
makes it particularly suitable for the study of the effect of 
baryonic physics on cosmic statistics. Such an intermediate 
box size, however, will be affected at both small and large 
scales by resolution and finite volume effects. 

At very small scales, counts-in-cells measurements are 
expected to suffer from shot noise, as the density field is sam- 
pled by a finite number of particles. This translates into both 
statistical variance and bias at small scales, if using naive 
statistical estimators for the moments of the density field. 
Assuming particles trace the density field as a local Poisson 
process, it can be shown, however, that factorial moments 
defin e d in Eq. [5] are unbiased estimators l|Szapudi fc Szalavl 
1 19931 : iBcrnard eau et aT]|2002f ). We thus do not expect our 
measurements to be affected by Poisson noise at small scales. 

On large scales, the results will be contaminated by 
cosmic variance and finite volume effects. In cosmological 
simulations, statistical quantities are usually computed by 
taking the spatial average - instead of ensemble average - 
of local quantities over the single simulated volume. This 
prescription is only appropriate for scales corresponding to 
wavenumbers k for which the simulation provides sufficient 
independent samples. For a box of a given size L, the sam- 
pling of large scales k with t~ approaching L suffers from the 
decreasing number of independent modes. The low number 
of modes at low wavenumbers introduces variance on large 



scale quantities, which is purely statistical in nature. As pre- 
sented in earlier, we have measured the cosmic variance of 
the whole MareNostrum box, and while conservative esti- 
mates for the errors range from 5%-30% depending on the 
statistic and estimation method, it is our understanding that 
cosmic variance does not fundamentally affect our result. As 
previously mentioned, we have minimized the effect of cos- 
mic errors on our conclusions by only considering ratios of 
statistical quantities from simulations run with the same set 
of random phases. 

We have shown that, although different halo profiles can 
describe variance amplification due to dissipative physics at 
small scale by merely modifying the concentration param- 
eter c(M), the third moment S3 introduces additional con- 
straints on the inner profiles which cannot be reproduced 
by changing c(M) alone. The distinctive slope of Ss(R) at 
small scales seems characteristic of a higher mass concen- 
tration towards the core than NFW. Unsurprisingly, profiles 
with a core or relatively weak central density peaks do not 
describe well the effective total matter distribution in the 
presence of baryons. 

Instead, we have found that using a superposition of 
a NFW profile and an exponential profile yields realistic 
variance amplification and S3 gain for reasonable values of 
the concentration parameters Co and b, disk mass fraction 
fd and disk scale Ad. One should note that the values of 
the best-fitting Ad and fd parameters are in good agreement 
with the expected physical properties of the galaxies of the 
simulated MareNostrum universe. The fd = 0.09 value is 
quite compatible with o bserved and predicted baryon disk 
mass fractions (see, e.g.. ISomerville et al.1120081 ). 

In this model, we chose not to introduce any mass or 
redshift dependence in fd and Ad- For the latter, this as- 
sumption is sup ported in part by the weak dependence on 
mass of rd/Rvii (|Somerville et al1 [2008). On the other hand, 
for fd, a proper model should account for the variation 
of M/L as a functio n of halo mass in real galaxies (see, 
e.g.. I Yang etalll2003l ). We postpone this more elaborate ap- 
proach to a future paper. 

The modified Co and b parameters of Eq. (23 for the 
NFW profile correspond to a more concentrated CDM com- 
ponent than in the pure DMO model. This is in accordance, 
both qualit a tively and quantitatively, with the results of 
iRudd et al.l <j2008l ) . It is interesting to note that the vari- 
ance boost caused by the NFW component is of the same 
order of magnitude as the adiabatic contraction effect visible 
on Fig. [2l albeit slightly weaker. This supports the idea that 
the composite halo profile concentrates a significant fraction 
of the halo's baryonic mass within the central disk, while the 
remaining halo gas essentially follows the NFW component 
which accounts for the cold dark matter. This last CDM 
component "feels" the presence of the hot gas through the 
process of adiabatic contraction. 

This suggests that both variance and skewness of the 
density field can be estimated at small scale within the 
framework of the halo model by using a composite halo pro- 
file. While the halo model is a valuable tool for the study 
of theoretical power spectra, the accuracy requirements of 
precision cosmology are arguably too stringent to consider 
directly fitting cosmological parameters to observed cosmic 
statistics. The halo model has merit however, as it allows 
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one to study the dependence of cosmic shear with baryonic 
features such as galaxy disk masses and sizes. 
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